#include "DXUT.h"
#include "ProjectedPolyhedron_d.h"
#include <stack>
using namespace std;

#define EPS 1.0e-3

PP::PP(void)
{
	nEquation = 0;
	nDegXYZ = NULL;
	coeff = NULL; 
}

PP::~PP(void)
{
	if (nDegXYZ)
	{
		delete []nDegXYZ;
		nDegXYZ = NULL;
	}

	if (coeff)
	{
		delete coeff;
		coeff = NULL;
	}
}

void PP::SetEquationCoeff(const int Equa, int* Deg, int maxDeg[3], double* weight)
{
	if (nDegXYZ)
	{
		delete[] nDegXYZ;
		nDegXYZ = NULL;
	}
	nDegXYZ = new int[3*Equa];
	for (int i=0; i<3*Equa; i++)
		nDegXYZ[i] = Deg[i];

	int ind[4];
	nEquation  = ind[0] = Equa;
	ind[1] = maxDeg[0]+1;
	ind[2] = maxDeg[1]+1;
	ind[3] = maxDeg[2]+1;
	nMaxDeg[0] = maxDeg[0];
	nMaxDeg[1] = maxDeg[1];
	nMaxDeg[2] = maxDeg[2];

	nMaxXYZDeg = max(maxDeg[0], maxDeg[1]);
	nMaxXYZDeg = max(nMaxXYZDeg, maxDeg[2]);

	if (coeff)
	{
		delete coeff;
		coeff = NULL;
	}
	coeff = new MultiArr<double>(4, ind, weight);
}

void PP::SolveEquations(int* RootNum, float* RootXYZ)
{
	double interval[3][2] = {{0.0f,1.0f}, {0.0f,1.0f}, {0.0f,1.0f}};
	MultiArr<double>* Coeff = new MultiArr<double>(*coeff);
	(*RootNum) = 0;

	//Iteration solution, if want test local interval, first ShrinkInterval.
	//ShrinkInterval(interval, Coeff);
	SolveEquationsIter(interval, Coeff, RootNum, RootXYZ);

	//No iteration solution!
	//SolveEquationsNoIter(interval, Coeff, RootNum, RootXYZ);

	delete Coeff;
	Coeff = NULL;
}

void PP::SolveEquationsNoIter(double Interval[3][2], MultiArr<double>* Coeff, int* RootNum, float* RootXYZ)
{
	typedef struct  
	{
		float Inter[3][2];
	} intervalXYZ;
	//intervalXYZ IntervalXYZ;
	stack<intervalXYZ> IntervalStack;
	stack<intervalXYZ>::value_type IntervalXYZ;

	double ResultLength[3], Scale[3];
	double AxisInterval[3][2] = {{0.0f,1.0f}, {0.0f,1.0f}, {0.0f,1.0f}};

	IntervalXYZ.Inter[0][0] = Interval[0][0]; IntervalXYZ.Inter[0][1] = Interval[0][1];
	IntervalXYZ.Inter[1][0] = Interval[1][0]; IntervalXYZ.Inter[1][1] = Interval[1][1];
	IntervalXYZ.Inter[2][0] = Interval[2][0]; IntervalXYZ.Inter[2][1] = Interval[2][1];
	IntervalStack.push(IntervalXYZ);

	while (!IntervalStack.empty() && (*RootNum)<10)
	{
		IntervalXYZ = IntervalStack.top();
		IntervalStack.pop();
		Interval[0][0] = IntervalXYZ.Inter[0][0]; Interval[0][1] = IntervalXYZ.Inter[0][1];
		Interval[1][0] = IntervalXYZ.Inter[1][0]; Interval[1][1] = IntervalXYZ.Inter[1][1];
		Interval[2][0] = IntervalXYZ.Inter[2][0]; Interval[2][1] = IntervalXYZ.Inter[2][1];

		ShrinkInterval(Interval, Coeff);

		GetProjectedIntevalIntersection(Coeff, Interval, AxisInterval); 

		if (AxisInterval[0][0]==1.0f || AxisInterval[1][0]==1.0f || AxisInterval[2][0]==1.0f)
		{
			continue;
		}

		for (int i=0; i<3; i++)
		{
			Scale[i] = AxisInterval[i][1]-AxisInterval[i][0];
			ResultLength[i] = Scale[i]*(Interval[i][1]-Interval[i][0]);
			AxisInterval[i][0] = (1-AxisInterval[i][0])*Interval[i][0] + AxisInterval[i][0]*Interval[i][1];
			AxisInterval[i][1] = (1-AxisInterval[i][1])*Interval[i][0] + AxisInterval[i][1]*Interval[i][1];
		}

		if (ResultLength[0]<EPS && ResultLength[1]<EPS && ResultLength[2]<EPS)
		{
			RootXYZ[3*(*RootNum)+0] = (AxisInterval[0][0]+AxisInterval[0][1])/2.0f;
			RootXYZ[3*(*RootNum)+1] = (AxisInterval[1][0]+AxisInterval[1][1])/2.0f;
			RootXYZ[3*(*RootNum)+2] = (AxisInterval[2][0]+AxisInterval[2][1])/2.0f;
			(*RootNum) = (*RootNum)+1;
			continue;
		}

		if ((Scale[0]>0.95&&ResultLength[0]>EPS) || (Scale[1]>0.95&&ResultLength[1]>EPS) || (Scale[2]>0.95&&ResultLength[2]>EPS) )
		{
			//intervalXYZ IntervalXYZRight;
			stack<intervalXYZ>::value_type IntervalXYZRight;

			//Subdivision
			bool NotSplit = true;
			for(int i=0; i<3; i++)
			{
				if(Scale[i]>0.95 && ResultLength[i]>EPS && NotSplit)
				{
					IntervalXYZ.Inter[i][0] = AxisInterval[i][0];
					IntervalXYZRight.Inter[i][1] = AxisInterval[i][1];
					IntervalXYZ.Inter[i][1] = IntervalXYZRight.Inter[i][0] = (IntervalXYZ.Inter[i][0]+IntervalXYZRight.Inter[i][1])/2;
					NotSplit = false;
				}
				else
				{
					IntervalXYZ.Inter[i][0] = IntervalXYZRight.Inter[i][0] = AxisInterval[i][0];
					IntervalXYZ.Inter[i][1] = IntervalXYZRight.Inter[i][1] = AxisInterval[i][1];
				}
			}

			IntervalStack.push(IntervalXYZ);
			IntervalStack.push(IntervalXYZRight);
		}
		else
		{
			for (int i=0; i<3; i++)
			{
				IntervalXYZ.Inter[i][0] = AxisInterval[i][0];
				IntervalXYZ.Inter[i][1] = AxisInterval[i][1];
			}
			IntervalStack.push(IntervalXYZ);
		}		
	}

	if (!IntervalStack.empty())
	{
		printf("Silhouette point number overflow!\n");
	}

}



void PP::SolveEquationsIter(double Interval[3][2], MultiArr<double>* Coeff, int* RootNum, float* RootXYZ)
{
	if ((*RootNum)>=10)
	{
		printf("Silhouette point number overflow!\n");
		return;
	}

	double ResultLength[3], Scale[3];
	double AxisInterval[3][2] = {{0.0f,1.0f}, {0.0f,1.0f}, {0.0f,1.0f}};

	GetProjectedIntevalIntersection(Coeff, Interval, AxisInterval); 

	if (AxisInterval[0][0]==1.0f || AxisInterval[1][0]==1.0f || AxisInterval[2][0]==1.0f)
	{
		return;
	}

	for (int i=0; i<3; i++)
	{
		Scale[i] = AxisInterval[i][1]-AxisInterval[i][0];
		ResultLength[i] = Scale[i]*(Interval[i][1]-Interval[i][0]);
		AxisInterval[i][0] = (1-AxisInterval[i][0])*Interval[i][0] + AxisInterval[i][0]*Interval[i][1];
		AxisInterval[i][1] = (1-AxisInterval[i][1])*Interval[i][0] + AxisInterval[i][1]*Interval[i][1];
	}

	if (ResultLength[0]<EPS && ResultLength[1]<EPS && ResultLength[2]<EPS)
	{
		RootXYZ[3*(*RootNum)+0] = (AxisInterval[0][0]+AxisInterval[0][1])/2.0f;
		RootXYZ[3*(*RootNum)+1] = (AxisInterval[1][0]+AxisInterval[1][1])/2.0f;
		RootXYZ[3*(*RootNum)+2] = (AxisInterval[2][0]+AxisInterval[2][1])/2.0f;
		(*RootNum) = (*RootNum)+1;
		return;
	}
	

	if ((Scale[0]>0.95&&ResultLength[0]>EPS) || (Scale[1]>0.95&&ResultLength[1]>EPS) || (Scale[2]>0.95&&ResultLength[2]>EPS) )
	{
		MultiArr<double>* CoeffRight = new MultiArr<double>(*Coeff);
		double IntervalRight[3][2];

		//Subdivision
		bool NotSplit = true;
		for(int i=0; i<3; i++)
		{
			if(Scale[i]>0.95 && ResultLength[i]>EPS && NotSplit)
			{
				Interval[i][0] = AxisInterval[i][0];
				IntervalRight[i][1] = AxisInterval[i][1];
				Interval[i][1] = IntervalRight[i][0] = (Interval[i][0]+IntervalRight[i][1])/2;
				NotSplit = false;
			}
			else
			{
				Interval[i][0] = IntervalRight[i][0] = AxisInterval[i][0];
				Interval[i][1] = IntervalRight[i][1] = AxisInterval[i][1];
			}
		}
		ShrinkInterval(IntervalRight, CoeffRight);	
		ShrinkInterval(Interval,      Coeff     );
		
		SolveEquationsIter(Interval,      Coeff,      RootNum, RootXYZ);
		SolveEquationsIter(IntervalRight, CoeffRight, RootNum, RootXYZ);
		delete CoeffRight;
		CoeffRight = NULL;
	}
	else
	{
		for (int i=0; i<3; i++)
		{
			Interval[i][0] = AxisInterval[i][0];
			Interval[i][1] = AxisInterval[i][1];
		}
		ShrinkInterval(Interval, Coeff);
		SolveEquationsIter(Interval, Coeff, RootNum, RootXYZ);
	}
}

void PP::ShrinkInterval(double subdivisionInterval[3][2], MultiArr<double>* weight)
{
	int otherXYZ[3][2] = {{1,2}, {0,2}, {0,1}};
	int ind[4], deg;	
	double *temp = new double[nMaxXYZDeg+1];
	double u;
	weight->SetData(coeff->GetData());
	for (int i=0; i<3; i++)
	{
		for (int j=0; j<nEquation; j++)
		{
			ind[0] = j;
			for (int m=0; m<=nDegXYZ[j*3+otherXYZ[i][0]]; m++)
			{
				ind[otherXYZ[i][0]+1] = m;
				for (int n=0; n<=nDegXYZ[j*3+otherXYZ[i][1]]; n++)
				{
					ind[otherXYZ[i][1]+1] = n;

					//Fetch the original weight;
					deg = nDegXYZ[j*3+i];
					for (int k=0; k<=deg; k++)
					{
						ind[i+1] = k;
						temp[k] = weight->FEI(ind);
					}

					//Subdivision [t0 t1] segment in [0, 1];
					//First use deCasteljaul fetch the right segment of t0;
					u = subdivisionInterval[i][0];
					for (int k=1; k<=deg; k++)
					{							
						for (int kk=0; kk<=deg-k; kk++)
						{
							temp[kk] = (1.0f-u)*temp[kk] + u*temp[kk+1];
						}
					}
					//Second, fetch the left (t1-t0)/(1-t0) segment of above;
					u = (subdivisionInterval[i][1]-u)/(1.0f-u);
					for (int k=1; k<=deg; k++)
					{
						for (int kk=k; kk<=deg; kk++)
						{
							temp[kk] = (1.0f-u)*temp[kk-1] + u*temp[kk];
						}
					}
					
					//Save the result;
					for (int k=0; k<=deg; k++)
					{ 
						ind[i+1] = k;
						weight->SEI(ind, temp[k]);
					}

				}//end of n
			}//end of m
		}//end of j
	}//end of i

	delete[] temp;
	temp = NULL;


}

void PP::GetProjectedIntevalIntersection(MultiArr<double>* Coeff, double Interval[3][2], double AxisInterval[3][2])
{
	int ind[4], deg;
	int otherXYZ[3][2] = {{1,2}, {0,2}, {0,1}};
	double maxCoeff, minCoeff, temp;	
	double* convex = new double[(nMaxXYZDeg+1)*2];
	double axisConvexIntersection[2] = {0.0f, 1.0f};

	AxisInterval[0][0] = 0.0; AxisInterval[0][1] = 1.0; 
	AxisInterval[1][0] = 0.0; AxisInterval[1][1] = 1.0; 
	AxisInterval[2][0] = 0.0; AxisInterval[2][1] = 1.0; 


	for (int i=0; i<3; i++)
	{
		axisConvexIntersection[0] = 0.0f;
		axisConvexIntersection[1] = 1.0f;
		for (int j=0; j<nEquation; j++)
		{
			ind[0] = j;
			deg = nDegXYZ[j*3+i];
			for (int k=0; k<=deg; k++)
			{
				ind[i+1] = k;
				maxCoeff = -FLT_MAX;
				minCoeff = FLT_MAX;
				for (int m=0; m<=nDegXYZ[j*3+otherXYZ[i][0]]; m++)
				{
					ind[otherXYZ[i][0]+1] = m;
					for (int n=0; n<=nDegXYZ[j*3+otherXYZ[i][1]]; n++)
					{
						ind[otherXYZ[i][1]+1] = n;	
						temp = Coeff->FEI(ind);
						maxCoeff = max(temp, maxCoeff);
						minCoeff = min(temp, minCoeff);
						//printf("%10.7f, %25.16e\n", temp, maxCoeff);
						//maxCoeff = temp>maxCoeff?temp:maxCoeff;
						//printf("%10.7f, %25.16e\n", temp, maxCoeff);
					}
				}
				convex[k]       = maxCoeff;
				convex[deg+1+k] = minCoeff;					
			}
			if (AxisIntersection(convex, axisConvexIntersection, deg))
			{
				AxisInterval[i][0] = max(axisConvexIntersection[0], AxisInterval[i][0]);
				AxisInterval[i][1] = min(axisConvexIntersection[1], AxisInterval[i][1]);
				if (AxisInterval[i][0]>AxisInterval[i][1])
				{
					AxisInterval[i][0] = 1.0f;
				}
			}

		}//end of j
	}//end of i

	delete[] convex;
	convex = NULL;
}

bool PP::AxisIntersection(double* convex, double intersect[2], int deg)
{
	//First Step: construct the convex using the max and min branch;
	int *convexPointIndex = new int[deg+1];
	int convexPoint = 1;	
	convexPointIndex[0] = 0;
	convexPointIndex[1] = 1;

	//max branch
	for (int i=2; i<=deg; i++)
	{
#if	UNIFORM_KNOT
		while (convex[i]-convex[convexPoint]>=convex[convexPoint]-convex[convexPoint-1] && convexPoint>=1)
#else
#endif
		{
			convexPoint--;
		}
		convex[++convexPoint] = convex[i];
		convexPointIndex[convexPoint] = i;
	}
	for (int i=1; i<=convexPoint; i++)
	{
		if (convex[i-1]<0.0f&&convex[i]>0.0f)
		{
#if	UNIFORM_KNOT
			intersect[0] = max(intersect[0], (convexPointIndex[i-1] + (convexPointIndex[i]-convexPointIndex[i-1])*(-convex[i-1])/(convex[i]-convex[i-1]))/deg);
#else
#endif
		}
		else if (convex[i-1]>0.0f&&convex[i]<0.0f)
		{
#if	UNIFORM_KNOT
			intersect[1] = min(intersect[1], (convexPointIndex[i-1] + (convexPointIndex[i]-convexPointIndex[i-1])*(-convex[i-1])/(convex[i]-convex[i-1]))/deg);
#else
#endif
		}
	}

	if (intersect[0]==0.0f && intersect[1]==1.0f && (convex[0]<0.0f))
	{
		delete[] convexPointIndex;
		intersect[0] = 1.0f;
		return false;
	}


	//min branch
	convexPoint = 1;	
	convexPointIndex[0] = 0;
	convexPointIndex[1] = 1;
	for (int i=0; i<=deg; i++)
	{
		convex[i] = convex[deg+1+i];
	}
	for (int i=2; i<=deg; i++)
	{
#if	UNIFORM_KNOT
		while (convex[i]-convex[convexPoint]<=convex[convexPoint]-convex[convexPoint-1] && convexPoint>=1)
#else
#endif
		{
			convexPoint--;
		}
		convex[++convexPoint] = convex[i];
		convexPointIndex[convexPoint] = i;
	}
	for (int i=1; i<=convexPoint; i++)
	{
		if (convex[i-1]<0.0f&&convex[i]>0.0f)
		{
#if	UNIFORM_KNOT
			intersect[1] = min(intersect[1], (convexPointIndex[i-1] + (convexPointIndex[i]-convexPointIndex[i-1])*(-convex[i-1])/(convex[i]-convex[i-1]))/deg);
#else
#endif
		}
		else if (convex[i-1]>0.0f&&convex[i]<0.0f)
		{
#if	UNIFORM_KNOT
			intersect[0] = max(intersect[0], (convexPointIndex[i-1] + (convexPointIndex[i]-convexPointIndex[i-1])*(-convex[i-1])/(convex[i]-convex[i-1]))/deg);
#else
#endif
		}
	}

	delete[] convexPointIndex;
	if (intersect[0]==0.0f && intersect[1]==1.0f && (convex[0]>0.0f))
	{
		intersect[0] = 1.0f;
		return false;
	}

	return true;
}
